********************************************************************************
* Modelling of representative dataset from Air Products
* Lei Zhang
* 11/11/2013
********************************************************************************

SETS
    i   "Set of Air Products' plants"       /P1*P5/
    ii  "Set of Competitor's plants"        /PC1*PC4/
    j   "Set of Demand points(markets)"     /D1*D10/
    k   "Set of products"                   /O1/
    t   "Set of time periods"               /T1*T5/;

ALIAS(t,tt);

*-------------------------------------------------------------------------------

SCALAR M    "Big M KKT"               /10000/;
SCALAR CMAX "Max amount of expansion" /1000/;
SCALAR CMIN "Min amount of expansion" /1/;

$include tabledata

* TABLE DE(j,k,t)        "Demand of product each market at each time period(100ton/90day)"
* TABLE DIS(i,j)         "Distance between Air Products' plants and markets(mile)"
* TABLE DISC(ii,j)       "Distance between competitor's plants and markets(mile)"
* TABLE PRICE(k,t)       "Price of products at different time periods(MM$/100ton)"
* TABLE FALPHA(i,k,t)    "Fixed cost(MM$)"
* TABLE FBETA(i,k,t)     "Variable cost of Air Products(MM$/100ton)"
* TABLE FGAMMA(i,k,t)    "Expansion cost of Air Products(MM$/100ton)"
* TABLE FGAMMAH(i,k,t)   "First 100ton higher expansion cost part of Air Products(MM$/100ton)"
* TABLE FT(i,k,t)        "Distribution cost of Air Products(MM$/100ton/mile)"
* TABLE FCT(ii,k,t)      "Distribution cost of competitor(MM$/100ton/mile)"

FGAMMA(i,k,t) = FGAMMA(i,k,t) / 10;
FGAMMAH(i,k,t) = FGAMMAH(i,k,t) / 10;

*-------------------------------------------------------------------------------

VARIABLES
    npv             "Net present value"
    mu1(j,k,t)      "kkt"; 

BINARY VARIABLES
    x(i,k,t)        "Selection of capacity Air Products' expansion"
    z(i,k,t)        "If first expansion happend (if c > 0), z=sum(w)"
    w(i,k,t)        "When first expansion happend"
*    z1(j,k,t)       "kkt"
    z2(i,k,t)       "kkt"
    z3(ii,k,t)      "kkt"
;

POSITIVE VARIABLES
    y(i,j,k,t)      "Amount of product that Air Products sells to market"
    yc(ii,j,k,t)    "Amount of product that competitor sells to market"
    c(i,k,t)        "Capacity of Air products facility"
    dc(i,k,t)       "Capacity expansion of Air products"
    
    
    mu2(i,k,t)     "kkt"
    mu3(ii,k,t)    "kkt";

    c.fx('P1','O1','T1') = 72;
    c.l('P1',k,t) = 72;
    c.fx('P2','O1','T1') = 225;
    c.l('P2',k,t) = 225;
    c.fx('P3','O1','T1') = 450;
    c.l('P3',k,t) = 450;
    c.fx('P4','O1','T1') = 0;
    c.l('P4',k,t) = 0;
    c.fx('P5','O1','T1') = 0;
    c.l('P5',k,t) = 0;

PARAMETER cc(ii,k);
    cc('PC1','O1') = 270;
    cc('PC2','O1') = 450;
    cc('PC3','O1') = 90;
    cc('PC4','O1') = 630;

    w.fx('P1',k,'T1') = 1;
    w.fx('P2',k,'T1') = 1;
    w.fx('P3',k,'T1') = 1;

*-------------------------------------------------------------------------------

EQUATIONS
    obj             "Objective function maximize Air Products NPV"
    epd(i,k,t)      "Investment decision in Air Products capacity expansion"
    upb(i,k,t)      "Upper bound of capacity expansion"
    lob(i,k,t)      "Lower bound of capacity expansion"
    tmd(j,k,t)      "Demand satisfaction for all markets"
    spd(i,k,t)      "Capacity and supply of Air Products"
    spdc(ii,k,t)    "Capacity and supply of competitor"
    b1(i,k,t)       "When expansion exist, first expansion exist"
    b2(i,k,t)       "relation between z and w"
    kkt11(i,j,k,t)  "kkt"
    kkt12(ii,j,k,t) "kkt"
*    kkt21(j,k,t)    "kkt"
*    kkt22(j,k,t)    "kkt"
    kkt31(i,k,t)    "kkt"
    kkt32(i,k,t)    "kkt"
    kkt41(ii,k,t)
    kkt42(ii,k,t)
;

* npv = total_income - fixed_cost - variable_cost - expansion_cost - transportation_cost
obj..               npv =e= SUM((i,j,k,t), PRICE(k,t)*y(i,j,k,t))
                        - SUM((i,k,t), FALPHA(i,k,t)*z(i,k,t))
                        - SUM((i,k,t), FBETA(i,k,t)*c(i,k,t))
                        - SUM((i,k,t), FGAMMAH(i,k,t)*w(i,k,t)+FGAMMA(i,k,t)*dc(i,k,t));

* c(i,k,t) = c(i,k,t-1) + dc(i,k,t) (when t < card(t)-1)
epd(i,k,t)$(ord(t)+1 le card(t)).. c(i,k,t+1) =e= c(i,k,t) + dc(i,k,t);
* dc(i,k,t) <= CMAX*x(i,k,t)
upb(i,k,t)..        dc(i,k,t) =l= CMAX*x(i,k,t);
* dc(i,k,t) >= CMIN*x(i,k,t)
lob(i,k,t)..        dc(i,k,t) =g= CMIN*x(i,k,t);
* sumi_y + sumi_yc = market_j's_demand
tmd(j,k,t)..        SUM(i, y(i,j,k,t)) + SUM(ii, yc(ii,j,k,t)) =e= DE(j,k,t);
* sumj_y < c(i,k,t)
spd(i,k,t)..        c(i,k,t) =g= SUM(j, y(i,j,k,t));
* sumj_yc < cc(ii,k,t)
spdc(ii,k,t)..      cc(ii,k) =g= SUM(j, yc(ii,j,k,t));

b1(i,k,t)..         x(i,k,t) =l= z(i,k,t);
b2(i,k,t)..         z(i,k,t) =e= SUM(tt$(ord(t) ge ord(tt)), w(i,k,tt));

kkt11(i,j,k,t)..    PRICE(k,t) + FT(i,k,t)*DIS(i,j) + mu1(j,k,t) + mu2(i,k,t) =l= 0.0001;
kkt12(ii,j,k,t)..   PRICE(k,t) + FCT(ii,k,t)*DISC(ii,j) + mu1(j,k,t) + mu3(ii,k,t) =l= 0.0001;
*kkt21(j,k,t)..      mu1(j,k,t) =l= M * (1-z1(j,k,t));
*kkt22(j,k,t)..      SUM(i, y(i,j,k,t)) + SUM(ii, yc(ii,j,k,t)) - DE(j,k,t) =g= -M * z1(j,k,t);
kkt31(i,k,t)..      mu2(i,k,t) =l= M * (1-z2(i,k,t));
kkt32(i,k,t)..      SUM(j, y(i,j,k,t)) - c(i,k,t) =g= -M * z2(i,k,t);
kkt41(ii,k,t)..     mu3(ii,k,t) =l= M * (1-z3(ii,k,t));
kkt42(ii,k,t)..     SUM(j, yc(ii,j,k,t)) - cc(ii,k) =g= -M * z3(ii,k,t);

*-------------------------------------------------------------------------------

MODEL COMPETITION /all/;

option optcr = 0;

SOLVE COMPETITION USING mip MAXIMIZING npv;
